// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                        avtL2NormBetweenCurvesQuery.C                      //
// ************************************************************************* //

#include <avtL2NormBetweenCurvesQuery.h>

#include <Utility.h>

#include <string>
#include <vector>

// ****************************************************************************
//  Method: avtL2NormBetweenCurvesQuery::avtL2NormBetweenCurvesQuery
//
//  Programmer:   Hank Childs
//  Creation:     October 2, 2003
//
// ****************************************************************************

avtL2NormBetweenCurvesQuery::avtL2NormBetweenCurvesQuery()
{
}


// ****************************************************************************
//  Method: avtL2NormBetweenCurvesQuery::~avtL2NormBetweenCurvesQuery
//
//  Programmer:   Hank Childs
//  Creation:     October 2, 2003
//
// ****************************************************************************

avtL2NormBetweenCurvesQuery::~avtL2NormBetweenCurvesQuery()
{
}


// ****************************************************************************
//  Method: avtL2NormBetweenCurvesQuery::CompareCurves
//
//  Purpose:
//      Computes the L2Norm number of the input curves.
//
//  Notes:
//      For two curves, c_a and c_b, let f(x) = |c_a(x) - c_b(x)|.
//      Then the L2-Norm is defined to be the integral of f(x)^2 over x.
//      
//      The curves are already broken up into line segments.  But those line
//      segments have differing endpoints across the curves.  (Meaning curve
//      1 may have a line segment from 0-0.435 in X, but curve 2's has a line
//      segment that goes from 0-0.546.  There are many more permutations of
//      this -- curve 2 may devote 3 segments to an interval in X that curve
//      1 covers in 1 line segment.)
//     
//      The easiest way to handle this problem is to create new curves that
//      have the same number of line segments and have the property that line
//      segment i for the new curves covers exactly the same interval in X.
//     
//      Once we have done that, we can start to consider what the contribution
//      of each line segment is to the L2-Norm.  Over line segment i, the
//      function f(x) is a line.  The line has the form f(x) = mx+b.  We also
//      know two points on the line -- f(x1) = |c1(x1) - c2(x1)| and
//      f(x2) = |c1(x2) - c2(x2)|.  Using these points we can calculate m and b.
//     
//      Once we know m and b, we can calculate what f(x)^2 is.  It is a 
//      parabola.  Finally, once we know f(x)^2, we can calculate the integral 
//      from x1 to x2.  
//      It is:
//      m*m*x1*x1*x1/3 + mb*x1*x1 + b*b*x1 - m*m*x2*x2*x2/3 + mb*x2*x2 + b*b*x2 
//     
//      So our gameplan is to break up the curves into line segments that run
//      over the same X-intervals, then calculate the function f(x), which is
//      a line, along that segment, and then determine that line segment's 
//      total contribution to the integral.  Once we have considered all of
//      the line segments, we will have calculated the whole integral.
//     
//  Programmer:   Hank Childs
//  Creation:     October 3, 2003
//
// ****************************************************************************

double
avtL2NormBetweenCurvesQuery::CompareCurves(int n1, const float *x1,
                     const float *y1, int n2, const float *x2, const float *y2)
{
    std::vector<float>  newCurve1Vals;
    std::vector<float>  newCurve2Vals;
    std::vector<float>  usedX;
    PutOnSameXIntervals(n1, x1, y1, n2, x2, y2, usedX, newCurve1Vals,
                        newCurve2Vals);

    int total_n_pts = (int)usedX.size();
    double integral = 0;
    for (int i = 0 ; i < total_n_pts-1 ; i++)
    {
        // pt1 and pt2 will determine f(x) for the range usedX[i]-[i+1]
        float pt1[2], pt2[2];
        pt1[0] = usedX[i];
        pt1[1] = newCurve1Vals[i] - newCurve2Vals[i];
        pt2[0] = usedX[i+1];
        pt2[1] = newCurve1Vals[i+1] - newCurve2Vals[i+1];

        // Calculate the line in f(x) = mx+b notation. 
        // Note that pt1[0] != pt2[0] because we threw out duplicated points.
        double rise = pt2[1] - pt1[1];
        double run  = pt2[0] - pt1[0];
        double m = rise / run;
        double b = pt1[1] - m*pt1[0];
        integral += 
                m*m*pt2[0]*pt2[0]*pt2[0]/3. + m*b*pt2[0]*pt2[0] + b*b*pt2[0]
              - m*m*pt1[0]*pt1[0]*pt1[0]/3. - m*b*pt1[0]*pt1[0] - b*b*pt1[0];
    }

    return sqrt(integral);
}


// ****************************************************************************
//  Method: avtL2NormBetweenCurvesQuery::CreateMessage
//
//  Purpose:
//      Creates a message for the L2Norm query.
//
//  Programmer: Hank Childs
//  Creation:   October 3, 2003
//
//  Modifications:
//
//    Cyrus Harrison, Tue Sep 18 13:45:35 PDT 2007
//    Added support for user settable floating point format string
//
//    Kathleen Biagas, Thu Feb 13 15:04:58 PST 2014
//    Add Xml results.
//
// ****************************************************************************

std::string
avtL2NormBetweenCurvesQuery::CreateMessage(double l2norm)
{
    MapNode result_node;
    result_node["L2Norm_between_curves"] = l2norm;
    queryAtts.SetXmlResult(result_node.ToXML());

    char msg[1024];
    std::string format = "The L2Norm between the two curves is " 
                    + queryAtts.GetFloatFormat() +".";
    snprintf(msg,1024,format.c_str(), l2norm);
    std::string m = msg;
    return m;
}


